#ifndef _BASIS_STATE_
#define _BASIS_STATE_
#include "../DF_Common/com_sparse_row.h"
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <string>
#include <vector>
using std::cout;
using std::endl;
using std::vector;

struct channel_TP {
  int num;
  int J_2;
  int T_2;
  int T_2z;
  int E_min;
  int E_max;
  int E_min_F;
  int E_min_I;
  int E_max_F;
  int E_max_I; // calculate in build_conf_T
  vector<int> indexE;
  vector<int> indexE_beg;
};

struct channel_Jacobi {
  int num;
  int J_2;
  int T_2;
  int E_min;
  int E_max;

  vector<int> indexE; // was marked in Jacobi.select() and Jacobi_NL.sort();
  vector<int> indexE_beg;
};

struct State_SP {
  // single particle state
  int n;
  int l;
  int E;
  int j_2;
  int tz_2; // pn:-1: proton ; 1: neutron tz:0
  int index;
};

struct State_SPM {
  // single particle state
  int n;
  int l;
  int E;
  int j_2;
  int m_2;
  int tz_2; // pn:-1: proton ; 1: neutron tz:0
  int index;
};

struct State_TP {
  // triple particle state
//public:
  int a;
  int b;
  int c;
  int Jab_2;
  int Tab_2;
  int Tab_2z; // pn
  int J_2;
  int T_2;
  int T_2z;
  int E;
  int Par;
  int F;         // n1, n1>0: exist in final_state
  int I;         // n2, n2>0: exist in initial_state; n1 == n2 come from same
                 // TP_allow_file
  int Ecut_Flag; // default = -1; changed in build_conf_Tz -> 1;
  int index;
//   State_TP() {
//     a = 0;
//     b = 0;
//     c = 0;
//     Jab_2 = 0;
//     Tab_2 = 0;
//     Tab_2z = 0; // pn
//     J_2 = 0;
//     T_2 = 0;
//     T_2z = 0;
//     E = 0;
//     Par = 0;
//     F = -1; // n1, n1>0: exist in final_state
//     I = -1;
//     Ecut_Flag = -1;
//     index = -1;
//
// #pragma acc enter data copyin(this)
// #pragma acc enter data copyin(a, b, c, Jab_2, Tab_2, Tab_2z, J_2, T_2, T_2z,   \
//                               E, Par, F, I, Ecut_Flag, index)
//   }
//   State_TP(int a_t, int b_t, int c_t, int Jab_2_t, int Tab_2_t, int Tab_2z_t,
//            int J_2_t, int T_2_t, int T_2z_t, int E_t, int Par_t, int F_t,
//            int I_t, int Ecut_Flag_t, int index_t) {
//     a = a_t;
//     b = b_t;
//     c = c_t;
//     Jab_2 = Jab_2_t;
//     Tab_2 = Tab_2_t;
//     Tab_2z= Tab_2z_t; // pn
//     J_2 = J_2_t;
//     T_2 = T_2_t;
//     T_2z = T_2z_t;
//     E = E_t;
//     Par = Par_t;
//     F = F_t;       // n1, n1>0: exist in final_state
//     I = I_t;       // n2, n2>0: exist in initial_state; n1 == n2 come from same
//                    // TP_allow_file
//     Ecut_Flag=Ecut_Flag_t; // default = -1; changed in build_conf_Tz -> 1;
//     index=index_t;
// #pragma acc enter data copyin(this)
// #pragma acc enter data copyin(a, b, c, Jab_2, Tab_2, Tab_2z, J_2, T_2, T_2z,   \
//                               E, Par, F, I, Ecut_Flag, index)
//   }
//   ~State_TP() {
// #pragma acc exit data delete (a, b, c, Jab_2, Tab_2, Tab_2z, J_2, T_2, T_2z,   \
//                               E, Par, F, I, Ecut_Flag, index)
//   }
//   void update_device() {
// #pragma acc update device(a, b, c, Jab_2, Tab_2, Tab_2z, J_2, T_2, T_2z, E,    \
//                           Par, F, I, Ecut_Flag, index)
//   }
//   void update_host() {
// #pragma acc update host(a, b, c, Jab_2, Tab_2, Tab_2z, J_2, T_2, T_2z, E, Par, \
//                         F, I, Ecut_Flag, index)
//   }
};

struct State_Conf {
  int F;    // index of final state
  int I;    // index of initial state
  int E_ex; // excited energy | E_{omega} + E_{hole} - E_{particle} | < Ecut
};
struct State_TPconf_E {
  // combination of E in each channel;
  int E_F;
  int E_I;
  friend bool operator==(const State_TPconf_E &a, const State_TPconf_E &b) {
    return ((a.E_F == b.E_F) && (a.E_I == b.E_I));
  }
};
struct State_TPM {
  // triple particle state
  int a_m;
  int b_m;
  int c_m;
  int Jab_2m;
  int Tab_2z; // pn
  int J_2m;
  int T_2z;
  int E;
  int index;
};

struct State_Jacobi {
  int n12;   // -- 0
  int n3;    // -- 1
  int l12;   // -- 2
  int l3;    // -- 1
  int S12_2; // -- 1
  int J12_2; // -- 1
  int j3_2;  // -- 1
  int T12_2; // -- 1
  int J_2;   // -- 1
  int T_2;   // -- 1
  int Par;
  int JT_channel;
  int E;
  int F;     // 1: exist in final_state
  int I;     // 1: exist in initial_state
  int index; // was marked in Jacobi.select()
  int flag;  // 1 for completeness group of JT-E channel; -1 for incompleteness
  // for Jacbobi_NL
  int N_0;
  int L_0;
  int E_NL;
  int Par_NL;
  int J2_NL;
  int JTNL_channel;
  int index_NL; // was marked in Jacobi_NL.sort()
  int F_NL;     // 1: exist in final_state
  int I_NL;     // 1: exist in final_state
};

class SingleP {
  // single particle state
public:
  vector<State_SP> state;
  vector<State_SP> state_Tz;   // pn
  vector<State_SPM> state_TzM; // pn
  vector<int> index_M2J;       // state_TzM -> state_Tz
  int n_max, l_max;
  void build_Tz(int n_max_t, int l_max_t);

  void build_Tz(int n_max_t, int l_max_t, int l_bad,
                int j2_bad); // without l_bad j2_bad : pn
  void build_TzM(int n_max_t, int l_max_t, int l_bad,
                 int j2_bad); // without l_bad j2_bad : pn
  // void build_T(int n_max_t, int l_max_t);
  void build_T(int n_max_t, int l_max_t, int l_bad,
               int j2_bad); // without l_bad j2_bad : T
  void build_T(std::string filename);
  void build_T(vector<vector<int>> sp_vec);
  void build_T_E(std::string filename);
  void build_T_E(int E_max);
  void build_Tz(std::string filename);
  void build_Tz(); // use build_T to build state_Tz
  void print();
  void print(int a_min, int a_max, int b_min, int b_max, int c_min, int c_max);
  void print(vector<vector<int>> TP_allow_vec);
  void print(vector<vector<int>> TP_allow_vec, std::ofstream &f_file);
  void print(std::ofstream &f_file);
  void print_DP(vector<vector<int>> DP_allow_vec);
  void print_DP_Tz(vector<vector<int>> DP_allow_vec);
  void print(std::string filename);
  void print_Tz();
  void print_Tz(std::string filename);
  void print_TzM(std::string filename);
  int FindSP(int n, int l, int j2, int tz_2);
  int FindSP_Tz(int n, int l, int j2, int tz_2);
};

struct State_DP {
  // double particle state
  int a;
  int b;
  int J_2;
  int T_2;
  int T_2z; // pn
  int E;
  int Par; // 0 or 1
  int index;
  int I;
  int F;
};

struct Channel_DP {
  int num;
  int J_2;
  int T_2;
  int T_2z;
  int E_min;
  int E_max;
  vector<int> indexE;
  vector<int> indexE_beg;

  int E_min_F;
  int E_min_I;
  int E_max_F;
  int E_max_I;
};

class DoubleP {
  // double particle state
public:
  vector<vector<State_DP>> state;    // state[channel][i]
  vector<vector<State_DP>> state_Tz; // state_Tz[channel][i]

  vector<Channel_DP> channel;
  vector<Channel_DP> channel_Tz;

  vector<vector<State_Conf>> DPconf,
      DPconf_Tz; // [ch][i].index_F & .index_I
  // the configuration of two TP state was
  // built in   void build_conf();
  vector<vector<vector<State_Conf>>> DPconf_k,
      DPconf_k_Tz; // [k][ch][i].index_F & .index_I; k is the file No.
                   // the configuration of two TP state was
                   // built in   void build_conf();
  vector<vector<vector<State_TPconf_E>>> DPconf_E,
      DPconf_E_Tz; // [k][ch][i].E_F & .E_I
                   // the combation of E in two TP state
                   // built in   void build_conf();
  vector<vector<State_TPconf_E>> DPconf_E_tot; // [k][i].E_F & .E_I
                                               // built in   void build_conf();
  // int lmax,j2max,nmax;
  int J2_max;
  int J2_max_DP;
  int Emax;
  int nmax;
  int j2max_sp;
  int all_size_T;                                   // calculate in build_T
  int all_size_Tz;                                  // calculate in build_Tz
  int Emax_Final, Emax_Init, Emin_Final, Emin_Init; // calculated in Select_FI_T
  int E_cut_conf; // | E_{omega} + E_{hole} - E_{particle} | <
                  // Ecut_configuration
  // DoubleP(){this->all_size_T = 0;}
  int E_cut_Luigi, E_model_SP_Luigi;
  int _A;
  DoubleP() : _A(0){};
  void build_Tz(SingleP &sp_t, int a_beg, int a_end, int b_beg,
                int b_end); // pn
  void build_Tz_Half(SingleP &sp_t, int a_beg, int a_end, int b_beg,
                     int b_end); // pp,nn: a<b pn: |pn>
  void sort_Tz();
  int FindC_Tz(int J2, int T_2z);
  int FindDP_Tz(int J2, int T_2z, int a, int b);
  int FindDP_Tz(int &ch_Tz, int &a, int &b);
  int FindDP(int &ch, int &a, int &b);

  void build_Constrain(SingleP &sp_t);
  void build_Constrain(SingleP &sp_t, vector<int> ab_vec, vector<int> E_cut_vec,
                       vector<int> Same_particle_vec,
                       vector<vector<vector<int>>> dp_all_vec);
  void build_Constrain_Luigi(SingleP &sp_t, int E_modle_SP);
  void Select_FI_Tz(SingleP &sp_t, vector<vector<vector<int>>> allow_FI_gather);
  void Select_FI_Tz_Luigi(SingleP &sp_t, int E_cut, int E_modle_SP);
  void Delet_Sort_Tz();
  void build_conf_Tz(SingleP &sp_t, vector<vector<int>> E_Sign_gather,
                     vector<int> E_cut_gather,
                     vector<int> Same_Particles_gather);
  void Delet_conf_E_Tz();

  void print_Tz(std::string filename);
  void print_channel_Tz(std::string filename);

  // void build_T(SingleP &sp_t, int a_beg, int a_end, int b_beg, int b_end,
  //              int c_beg, int c_end); // T
  // void Select_FI_T(SingleP &sp_t, vector<vector<vector<int>>>
  // allow_FI_gather);
  // void build_conf_T(SingleP &sp_t, vector<vector<int>> E_Sign_gather,
  //                   vector<int> E_cut_gather,
  //                   vector<int> Same_Particles_gather);
  // void Delet_Sort_T();
  //
  // void Delet_conf_E_T();

  void print();
  void print(int channel_num);
  void print(std::string filename);

  void print_channel();
  void print_channel(std::string filename);
  void print_TPconf_info(std::string filename1, std::string filename2,
                         std::string filename3);
  void print_TPconf_info_Tz(std::string filename1, std::string filename2);
  void sort();
  int FindC(int J2, int T2);

  int FI_compare(int F, int I, int k);
  int FI_compare(int F, int I);

  int Tri(int &j2a, int &j2b, int &j2c);

private:
  int phase(int n) { return n % 2 ? -1 : 1; }
};

class TripleP {
  // triple particle state
public:
  vector<vector<State_TP>> state;    // state[channel][i]
  vector<vector<State_TP>> state_Tz; // state_Tz[channel][i]
  vector<State_TPM> state_TzM;       // state_TzM[i]
  vector<channel_TP> channel;
  vector<channel_TP> channel_Tz;
  vector<vector<int>> JT2channel; // calculate in build_JT();
                                  // JT2channel[0:1][1:J2max], 0: T=0.5; 1:T=1.5
  vector<vector<State_Conf>> TPconf,
      TPconf_Tz; // [ch][i].index_F & .index_I
  SMatrix_Ch TPconf_S, TPconf_Tz_S;
  // the configuration of two TP state was
  // built in   void build_conf();
  vector<vector<vector<State_Conf>>> TPconf_k,
      TPconf_k_Tz; // [k][ch][i].index_F & .index_I; k is the file No.
                   // the configuration of two TP state was
                   // built in   void build_conf();
  vector<vector<vector<State_TPconf_E>>> TPconf_E,
      TPconf_E_Tz; // [k][ch][i].E_F & .E_I
                   // the combation of E in two TP state
                   // built in   void build_conf();
  vector<vector<State_TPconf_E>> TPconf_E_tot; // [k][i].E_F & .E_I
                                               // built in   void build_conf();
  vector<vector<int>> E_Jac_conf;              //
  // vector<vector<int>> E_Jac_conf;
  // ----- for select build_T ------
  vector<vector<vector<int>>> TP_allow_gather_tp;

  // ------    ------   ------
  // int lmax,j2max,nmax;
  int J2_max;
  int J2_max_TP;
  int Emax;
  int nmax;
  int j2max_sp;
  int all_size_T;                                   // calculate in build_T
  int all_size_Tz;                                  // calculate in build_Tz
  int Emax_Final, Emax_Init, Emin_Final, Emin_Init; // calculated in Select_FI_T
  int E_cut_conf; // | E_{omega} + E_{hole} - E_{particle} | <
  // Ecut_configuration

  void init_Constrain(vector<vector<vector<int>>> TP_allow_gather);
  int flag_constrain(SingleP &sp_t, int a, int b, int c, int &Flag_F,
                     int &Flag_I);
  int flag_constrain_Tz(SingleP &sp_t, int a_Tz, int b_Tz, int c_Tz,
                        int &Flag_F, int &Flag_I) {
    int a = std::floor(a_Tz / 2);
    int b = std::floor(b_Tz / 2);
    int c = std::floor(c_Tz / 2);
    int res = flag_constrain(sp_t, a, b, c, Flag_F, Flag_I);
    return 0;
  }

  void build_Tz(SingleP &sp_t);             // pn
  void build_Tz(SingleP &sp_t, int Emax_t); // pn

  void build_Tz(SingleP &sp_t, int a_beg, int a_end, int b_beg, int b_end,
                int c_beg, int c_end);                // pn
  void build_Tz_Constrain(SingleP &sp_t, int Emax_t); // pn
  void build_Tz_Constrain(SingleP &sp_t, int a_beg, int a_end, int b_beg,
                          int b_end, int c_beg, int c_end); // pn
  void build_TzM(SingleP &SP, int Emax_t);                  // pn
  void build_T(SingleP &sp_t, int Emax_t);                  // T
  void build_T(SingleP &sp_t, int a_beg, int a_end, int b_beg, int b_end,
               int c_beg, int c_end);                // T
  void build_T_Constrain(SingleP &sp_t, int Emax_t); // T
  void build_T_Constrain(SingleP &sp_t, int a_beg, int a_end, int b_beg,
                         int b_end, int c_beg, int c_end); // T
  void build_TTz(SingleP &sp_t, int a_beg, int a_end, int b_beg, int b_end,
                 int c_beg, int c_end); // T

  void Select_FI_T(SingleP &sp_t, vector<vector<vector<int>>> allow_FI_gather);
  void Select_FI_Tz(vector<vector<vector<int>>> allow_FI_gather);
  void build_conf_T(SingleP &sp_t, vector<vector<int>> E_Sign_gather,
                    vector<int> E_cut_gather,
                    vector<int> Same_Particles_gather);
  void build_conf_T_S(SingleP &sp_t, vector<vector<int>> E_Sign_gather,
                      vector<int> E_cut_gather,
                      vector<int> Same_Particles_gather);
  void build_conf_Tz(SingleP &sp_t, vector<vector<int>> E_Sign_gather,
                     vector<int> E_cut_gather,
                     vector<int> Same_Particles_gather);
  void build_conf_Tz_S(SingleP &sp_t, vector<vector<int>> E_Sign_gather,
                       vector<int> E_cut_gather,
                       vector<int> Same_Particles_gather);
  void Delet_Sort_T();
  void Delet_Sort_Tz();
  void Delet_conf_T(SingleP &sp_t, vector<vector<int>> E_Sign_gather,
                    vector<int> E_cut_gather,
                    vector<int> Same_Particles_gather);
  void Delet_conf_Tz(SingleP &sp_t, vector<vector<int>> E_Sign_gather,
                     vector<int> E_cut_gather,
                     vector<int> Same_Particles_gather);

  void find_E_Jac_conf(vector<vector<int>> &E_cut_conf);

  void print();
  void print(int channel_num);
  void print(std::string filename);
  void print_Tz(std::string filename);
  void print_TzM(std::string filename);
  void print_channel();
  void print_channel(std::string filename);
  void print_channel_Tz(std::string filename);
  void print_TPconf_info(std::string filename1, std::string filename2,
                         std::string filename3);
  void print_TPconf_info_Tz(std::string filename1, std::string filename2);
  void sort();
  void sort_Tz();
  int FindC(int J2, int T2);
  int FindC_not_exit(int J2, int T2);
  int FindC_Tz(int J2, int T_2z);
  int FindTP_Tz(int J2, int T_2z, int J2ab, int a, int b, int c);
  int FindTP_not_exit(int ch, int J2ab, int T2ab, int a, int b, int c);
  int TzDiff(State_TP &A,
             State_TP &B); // if the only difference is Tz : 1; else -1;
  int FI_compare(int F, int I, int k);
  int FI_compare(int F, int I);

private:
  int Tri(int &j2a, int &j2b, int &j2c);

  int phase(int n) { return n % 2 ? -1 : 1; }
};

class Jacobi {
public:
  vector<vector<State_Jacobi>> state; // state[channel][i]
  vector<channel_Jacobi> channel;
  vector<vector<int>> JT2channel; // calculate in build_JT();
                                  // JT2channel[0:1][1:J2max], 0: T=0.5; 1:T=1.5
  // int lmax,j2max,nmax;
  int nmax, n12max, n3max; // calculate in sort()
  int lmax, l12max, l3max; // calculate in sort()
  int J12_2max, j3_2max;   // calculate in sort()
  int nmin, n12min, n3min; // calculate in sort()
  int lmin, l12min, l3min; // calculate in sort()
  int J12_2min, j3_2min;   // calculate in sort()
  int all_size;            // calculate in sort()
  int l12max_F, l3max_F, l12min_F, l3min_F;
  int l12max_I, l3max_I, l12min_I, l3min_I;
  vector<vector<int>>
      Jaco_l123_info; // [F/I][0-3] : l12min,l12max, l3min, l3max;
  int Emax_F, Emin_F, Emin_I, Emax_I;
  vector<int> Jaco_E_info; // E_F_min, E_F_max, E_I_min, E_I_max

  int JT_num;
  int J2max;
  int Emax;
  vector<vector<State_Conf>> Jacconf; // [ch][i].index_F & .index_I

  void build_JT(int n12max_t, int n3max_t, int l12max_t, int l3max_t);
  void build_JT(int Emax_t, int l12max_t, int l3max_t);
  void build_JT(int Emax_t);
  void build_JT(int Emax_t, int J2max_t);
  void sort();
  void print();
  void print(int channel_num);
  void print(std::string filename);
  void print_channel();
  void print_channel(std::string filename);
  void print_Jacconf(std::string filename);
  void select();
  void select_FI_T(vector<vector<State_TPconf_E>> TP_conf_E_tot);
  void build_conf_T(vector<vector<State_TPconf_E>> TP_conf_E_tot);
  void Delet_Sort_T(TripleP &TP);
  void Delete_unused_Jaco_channel(vector<int> Jaco_used_Flag);
  void check();

  int FindC(int J2, int T2);

protected:
  // int num;
  int Tri(int &j2a, int &j2b, int &j2c);
  int Tri_state(State_Jacobi &state);
  int phase(int n) { return n % 2 ? -1 : 1; }
};

class Jacobi_NL {
  // Jacobi with the quantum number (N0 L0) of C.M.
public:
  vector<vector<State_Jacobi>> state_NL; // state[channel][i]
  vector<channel_Jacobi> channel;
  vector<vector<int>> JT2channel; // calculate in build_JT();
                                  // JT2channel[0:1][1:J2max], 0: T=0.5; 1:T=1.5
  // int lmax,j2max,nmax;
  int N0_max, L0_max;
  int Emax;
  int J2abc_max;
  int all_size; // calculate in sort()
  int channel_size_T1;
  // channel_size_T1 calculate in sort(), used in FindC
  // the number of channel when T2=1, T2=3 : size - size_T1
  void build_NL(Jacobi &jaco, int N0_max_t, int L0_max_t);
  void build_NL(Jacobi &jaco, TripleP &TP);
  void sort();
  void print();
  void print(int channel_num);
  void print(std::string filename);
  void print_short(std::string filename);
  void print_channel();
  void print_channel(std::string filename);
  int FindC(int J2, int T2);
  void Delet_Sort_T();
  void select_FI_T(TripleP &TP);
  void Delete_unused_Jaco_channel(Jacobi &jaco, vector<int> &Jaco_used_Flag);

private:
  //
  int Tri(int &j2a, int &j2b, int &j2c);
  int phase(int n) { return n % 2 ? -1 : 1; }
};

struct TP_Toku {
  int channel;
  int T_2;
  int J_2;
  int index;
};

#endif
